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Modeling of random bimodal structures of 
composites (application to solid propellants): 
II. Estimation of effective elastic moduli 

V.A. Buryachenko 1 2) Q 

&i Abstract: We consider a linearly elastic composite medium, which consists of a homogeneous matrix 

containing a statistically homogeneous set of multimodal spherical inclusions modeling the morphology 



in 



of heterogeneous solid propellants (HSP). Estimates of effective elastic moduli are performed using the 
multiparticlc effective field method (MEFM) directly taking into account the interaction of different 
inclusions. Because of this, the effective elastic moduli of the HSP evaluated by the MEFM are sensitive 
to both the relative size of the inclusions (i.e., their multimodal nature) and the radial distribution 
functions (RDFs) estimated from experimental data, as well as from the ensembles generated by the 
method proposed. Moreover, the detected increased stress concentrator factors at the larger particles in 
comparison with smaller particles in bimodal structures is critical for any nonlinear localized phenomena 
for HSPs such as onset of yielding, failure initiation, damage accumulation, ignition, and detonation. 
Keywords: Microstructurcs, random structures, inhomogeneous material, effective clastic moduli 

1. Introduction 

Heterogeneous solid propellants (HSPs) commonly used in aerospace propulsion are likely to play 
an important role for as long as rockets are built. The typical HSPs consist of multimodal spherical 
particles of solid oxidizer (ammonium pcrchlorate, AP) of order 1-100 /im diam embedded in a rubbery 
fuel binder such as hydroxi-terminated-polybutadiene (HTPB) or polybutadiene acrylonitrile (PBAN) 
(see, e.g. Sutton and Biblarz, 2003). Designers of rockets are concerned with a number of propellant- 
related issues, including the burning rate, the thermal and mechanical properties of the propellants, and, 
for metallized propellants, the behavior of the metal particles at the surface, including agglomeration (see 
Massa, Jackson and Short, 2003; Matous and Geubelle, 2006; Matous et al, 2007; Maggi et al., 2008; Xu, 
Aravas and Sofronis, 2008). 

Both numerical and experimental investigations of random structures are intended for prediction of 
overall mechanical and physical properties of these media. Suspensions probably is one of the most exten- 
sively studied area of highly packed particulate microhetcrogencous media. Numerous authors minimize 
the non-hydrodynamic factors considering the low shear rate Newtonian limit that is most interesting 
from the point of view of analogy with composite materials. In simplest cases, one assumes that at the 
plotting of the relative effective viscosity r]*/r]^ of a suspension versus the adjustable volume fraction 
c/c m , viscosity versus concentration plots collapse onto one curve (see for references Stickel and Powell, 
2005; Arefinia and Shojaei, 2006): rf /rf°) = [1 - /(c«, . . . , c (™))]-^ c ™. He and Ekere (2001) (see also 
Chong, Christiansen and Baer, 1971; Shapiro and Probstein, 1992; Patlazhan, 1993; D'Haene and Mewis, 
1994) proposed semi- analytical viscosity model of noncolloidal bimodal suspension where c m depends 
upon both the size ratio A of large and small particles and the volume fraction £ of large particles. In a 
limiting case A — > oo, the mentioned models are reduced to a fruitful alternative by Farris (1968) (see 
also Probstein, Sengun and Tseng, 1994; Greenwood, Luckham and Gregory, 1998; Zaman and Moudgil, 
1998) based on purely geometric arguments. Assuming that the suspension of the fine particles can be 
treated as the suspending media for the large particles, Farris (1968) simulated the full bimodal mix- 
ture of particles by mapping it onto an effective one-modal system which can be simulated more easily. 
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In this model, the fine and coarse fractions are assumed to behave independently of each other and 
the viscosity of the mixture of smaller fractions is the medium viscosity for the next largest fraction: 
7y*/jy(°) = /(c (1) ) • /(c (2) ) • . . . • /(c (n) ) (a (1) <C a (2) < . . . < a (n) ). A similar approach (under the name 
a two-step homogenization technique) is well known (even if the name Farris (1968) is not mentioned) 
in the engineering material science of both the composites with the large size ratio A of fillers and ag- 
glomerated composite materials (for references see, e.g., Buryachenko, 2007a; Enikolopyan, et al, 1990; 
Clements and Mas, 2004). 

Composites with a wide range of bimodal inclusion concentrations are less well understood and 
call for further investigation. So, a standard so-called mult i- particle unit cell method of computational 
micromcchanics is based on pcriodization of generated random media (see for references Gusev, Lusti and 
Hine, 2002; Lusti and Gusev, 2004; Bohm, Han and Eckschlagcr, 2006; Buryachenko, 2007a) reducing to 
the modeling of periodic boundary conditions at the meso cell containing just a few tens of inclusions 
that eliminates both the boundary and size effects in a windowing approach (see, e.g., Liu et al., 2009). 
The generation process is based upon Monte Carlo simulation method wherein the particles are generated 
from a certain aggregate size distribution and then randomly placed into the meso cell in such a way 
that there is no intersection between the particles. Then, the generated mesostructure models are to be 
meshed for computational testing of the composite with the finite element method (FEM) in order to 
estimate of the macroscopic response [the use of the boundary element method (Chen and Papathanasiou, 
2004) and multipole method (Hatch and Davis, 2006; Kushch, Sevostianov and Mishnaevsky, 2008) is 
also known]. In so doing, the pack generations are based on either the hard core model simulation 
(Segurado and Llorca, 2002; Wisslcr et al., 2003; Hafner, et al, 2006; Wriggers and Moftah, 2006; Kushch, 
Sevostianov, and Mishnaevsky, 2008; Tawerghi and Yi, 2009) or the collective rearrangement model 
(CRM) destined for creation of the close packing (Stroeven, Askes, Sluys, 2004; Annapragada, Sun and 
Garimella, 2007; Matous, et al., 2007). Disadvantages of these random generation methods for composites 
with a wide range of particle concentration were considered by Buryachenko et al. (2003). The advances in 
imaging techniques [e.g., computed tomography (CT), scanning electron microscopy (SEM), and magnetic 
resonance imaging (MM)] require powerful computational methods for numerical analysis of CMs. To 
simplify the technology of finite clement mesh generation for particle reinforced material (e.g., HSP), 
enrichment techniques is used to account for the material interfaces in the framework of extended FEM 
(XFEM, see for references Fries and Bclytschko, 2010; Hiriyur, Waisman, and Dcodatis, 2011), first 
introduced by Bclytschko and Black (1999) as a solution to the rcmcshing issue for crack propagation. 
The geometry of material distribution is described by level set function, which allows one to model the 
internal boundaries of the microstructure without the adaptation of the mesh (Du, Ying and Jiang, 
2010; Legrain et al., 2011). Some alternative approaches based on a real structure can be done mainly 
in two different ways. The first one can directly import the experimentally obtained real structure 
into FE software (see for references Ghosh, 2011; Legrain et al., 2011). The second approach based on 
the concept of Statistically Equivalent Periodic Unit Cell (SEPUC) generates the model structures that 
have similar statistical functions as the experimentally obtained ones with a subsequent substitution of 
found descriptor into analytical micromcchanics models (see, e.g., Povirk, 1995; Matous, Leps, Zeman 
and Scjnoha, 2000; Zeman and Sejnoha, 2001; Borbely, Kenesei and Biermann, 2006; Lee, Gillman and 
Matous, 2011). 

The modeling of composites with bimodal spherical inclusions by the methods of analytical micromc- 
chanics (see for references Nemat-Nasser and Hori, 1993; Torquato, 2002; Milton, 2002; Buryachenko, 
2007a) has been developed less and necessitates some additional consideration. It is explained by a fact 
that the most popular methods of analytical micromechanics are essentially one-particle ones. First of all, 
this will be true for the variants of the effective medium method by Kroner (1958) and by Hill (1965) and 
the mean field method (Mori and Tanaka, 1973; Benveniste, 1987). Recently a new method has become 
known, namely the multiparticle effective field method (MEFM) was put forward and developed by one 
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of the authors (see for references Buryachenko, 2007a). The MEFM is based on the theory of functions of 
random variables and the Green functions. MEFM does not make use of a number of hypotheses which 
form the basis of the traditional one-particle methods. It is obvious that the classical methods mentioned 
above do not allow for direct binary interaction of inclusions, and, because of this, these methods are 
invariant to the size distribution of inclusions. Contrastingly, the MEFM has qualitative benefits following 
immediately from the consideration of multiparticle interactions and manifesting by a sensitivity of the 
MEFM's estimations to the size distribution of heterogeneities. 

So called discrete structural models alternatively taking binary interactions of inclusions into account 
are worthy of notice (see, e.g., Zgaevskii, 1977; Garishin, 1997; Garishin and Moshcv, 2002; Yanovsky and 
Zgaevskii, 2004). These methods were initiated by Chen and Acrivos (1978) who showed that the elastic 
energy stored in the matrix layer between the spheres is many times greater than that in the rest of 
the surrounding matrix. This matrix layer can be considered as the elastic rod or spring transferring the 
interaction force between the spheres that makes it possible to use the principle of physical discretization 
developed by Absi and Prager (1975). A composite material is modeled as an elastic network which seems 
to be much simpler for mathematical treatment than the original continuum mechanics representation. 
The fundamental limitations of the method are a simple composite structure and the growing of an 
approximate error of a spring model with both the growing of inter-particle distances and decreasing of 
clastic mismatch of constitutives. 

The paper is organized as follow. In section 2 we present the basic field equations of linear elasticity, 
notations, and statistical description of the composite microstructurc. In Section 3 we recall the basic 
concepts defining some classical method of micromechanics and present explicit formula for both effective 
clastic moduli and stress concentrator factor for composites with bimodal distribution of spherical parti- 
cles and the radial distribution functions (RDFs) evaluated in an accompanying paper by Buryachenko, 
Jackson and Amadio (2012), henceforth referred to as (I). In Section 4 it is detected that the effective 
elastic moduli of the HSP evaluated by the MEFM are very sensitive (that is confirmed by comparison 
with available experimental data) to both the relative size of inclusions (i.e. their multimodal nature) 
and the RDFs. Moreover, one detects increased stress concentrator factors at the big particles in com- 
parison with one for the small particles in bimodal structures that is critical for any nonlinear localized 
phenomena such as onset of yielding, failure initiation, damage accumulation, and others. A wide class 
of prospective problems for the modeling of HSPs is considered in Section 5. 

2. Preliminaries 

Let a linear elastic infinite body R d contains an open bounded domain w C R d (window of obser- 
vation) with a boundary T and with an indicator function W and space dimensionality d (d = 2 and 
d = 3 for 2-D and 3-D problems, respectively). The domain w contains a homogeneous matrix 
and a random finite set X = (vi) (i = 1, . . . ,N(w)) of inclusions Vi with centers x, and with charac- 
teristic functions Vi(y) equals one if y e Vi and zero otherwise and bounded by the spherical surfaces 



T^ 1 ' = {y : |y — x,| = a' 1 '} and I^ 2 ' = {y : |y — Xj| = a^} of two radii a*- 1 ) and a' 2 ', respectively; 
«(*) = Uvf } (1 = 1,2,..., fc = l,2). 



The local strain tensor e is related to the displacements u via the linearized strain-displacement 
equation 



Here (g> denotes tensor product, and (.) T denotes matrix transposition. The stress tensor a satisfies the 
equilibrium equation (no body forces acting): 



£ 




(1) 
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Stresses and strains are related to each other via the constitutive equations 

<t(x) = L(x)e(x) or e(x) = M(x)tr(x). (3) 

L(x) and M(x) = L(x)^ 1 are the known phase stiffness and compliance fourth-order tensors, and the 
common notations for scalar products have been employed: Le = L^kiSki- In particular, for isotropic 
constituents the local stiffness tensor L(x) is given in terms of the local bulk modulus fc(x) and the local 
shear modulus /i(x): 

L(x) = (dfc(x), 2/i(x)) = dfc(x)Ni + 2 M (x)N 2 , /3(x) = /3 (x) 8, (4) 
where in component form 

N\\ijki — -^SijS k i, N 2 \ijki = -^{SikSji + $uSjk) — -^ijSki (5) 

and Sij is the Kronecker delta; i,j, k, I = 1, 2, . . . , d; d = 2, 3. 

All tensors f (f = L,M) of material properties are decomposed as f = f c + fi(x). The subscript 1 
denotes a jump of the corresponding quantity (e.g. of the material tensor), f is assumed to be constant 
in the matrix tA°) = w\v and inside the inclusions: f(x) = f(°) for x e v^°\ and f(x) = f(°) + f[ for 
x 6 y( k \ Here and in the following the upper index V s ' numbers the components and the lower index 
i numbers the individual inclusions; v = Uv^ = Uvi, T^(x) = = S^( x )? an d W fc -*(x) is a 

characteristic function of v^ k \ equal to 1 if x G and otherwise, (k = 1,2; i = 1, 2, . . .). 

We assume that the phases are perfectly bonded, so that the displacements and the traction com- 
ponents of the stresses are continuous across the interphase boundaries. For the mesodomain x £ w we 
take either uniform displacement or traction boundary conditions 

u(x) = £°x, (6) 
cr(x)n(x) = cr n(x) = t(x), (7) 

respectively, where t(x) is the traction vector at the external boundary dw, n is its unit outward normal, 
and e° and <x° are the mesoscopic strain and stress tensors, i.e. a given constant symmetric tensor. 

All the random quantities under discussion are statistically homogeneous, and, hence, the ensemble 
averaging could be replaced by volume averaging 

((.)) =W- 1 J(.)W(x)d X , <(.)>« = [^ fc )]- 1 |(.)^W(x)dx, (8) 

where the bar appearing above the region represents its measure, e.g. v = mes v. The average over 
component agrees with the ensemble average over an individual inclusion Vi G (i = 1, 2, . . .) : 

((•))« = ii ))^ ■ The notation ((.))i(x) at x 6 Vi C v k means the average over an ensemble realization 
of surrounding inclusions (but not over the volume Vi of a particular inclusion, in contrast to ((.))i). 
For the description of the random statistically homogeneous structure of a composite material, let us 
introduce a conditional probability density f{vj, [; Vi, x^), which is a probability density to find the 
j-th inclusion with the center Xj in the domain Vj with fixed inclusion Vi with the centers x^, and 
x 3 ^ Xi. Of course, (p(vj , Xj | ; Uj, x,) for values of Xj lying inside the "excluded volumes" (called also the 
correlation hole) Ui>^, where D Vi with characteristic functions Vo m (since inclusions cannot overlap), 
and (/?(wj, Xj |; «i, x.;) — > ip(vj,x.j) at |xj — x.;| — > oo (since no long-range order is assumed); it is assumed 
that = Vj are the circles with the radii 2a; j, i = 1,2...). Only if the pair distribution function 



V.A. Buryachenko 



5 



gkmtx-j — Xj) = <p(vj, Xj|; Uj, Xi)/n^ m ^ (u^ G u^, G t/ m ', fc,m = 1,2) depends on |x« — Xj| it is called 
the radial distribution function. <p(yi,Xi) is a number density nS^ of component 9 Vi and c^-* is 
the concentration, i.e. volume fraction, of the component v^ k >: & k ' = (V^) = ViW- k > (fe = 1, 2; i = 
1,2,...), =l-(V),c=(V). 



3 Effective elastic moduli and stress concentrator factors 

We will slightly modify the basic assumptions and the final formula of the multiparticlc effective field 
method (MEFM) for estimation of effective elastic moduli of composites with polydispersc particles. This 
case was not considered by Buryachenko (2007a) where a detailed discussion and numerous references for 
this and related methods can be found. 

Substituting (3) and (1) into the equilibrium equation (2), we obtain a differential equation with 
respect to the displacement u which may be reduced to a symmetrized integral form for the stresses 

<t(x) = (*) + J r(x - y)fo(y) - (r,)]dy, (9) 

where the tensor rj{y) = Mi(y)er(y) is called the strain polarization tensor and is simply a notations! 
convenience (see Willis, 1981). The integral operator kernel 

T(x - y) = -1/°) [K(x - y) + U(x - y)l/ )] , (10) 

is the even homogeneous generalized functions of the order —d defined by the second derivative of the 
Green tensor G: 

U m (x) = [V^VAfcW]^), (11) 

where the notation indicates symmetrization on (ij) and (kl), I is a unit fourth-order tensor, and G is 
the infinite-homogeneous-body Green's function of the Navier equation with an elastic modulus tensor 
L(°), 

V [V ® G(x) + (V ® G(x)) T ] /2} = -W(x), (12) 

vanishing at infinity (|x| — > oo), 6(x) is the Dirac delta function, 6 is the unit second order tensor. The 
equation (12) is valid for the domain x G w containing statistically large number of inclusions except 
when the distance of x from the boundary dw of the body w is of order of a correlation length of the 
microstructure (see for details Willis, 1981; Buryachenko, 2007a). 

Both the effective compliance M* and the effective moduli L* = (M*)" 1 appearing in the overall 
constitutive relation 

(e) = M», (a) = L*(e) (13) 

can be defined by the general relation 

M* = (MB*), (14) 
where B* = B*(x) is a local stress concentration tensor obtained under pure mechanical loading 

< T(x)=B*(x)(a). (15) 



After conditional statistical averaging Eq. (9) turns into an infinite system of integral equations. 
In order to close and approximately solve this system we now apply the MEFM hypotheses (see details 
Byryachenko, 2007a): 
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HI) Each inclusion Vi has a spherical form and is located in the effective field CT,(y) which is 
homogeneous over the inclusion Vi and 

Wi(y)=Wi=^ (16) 

for any y £ Vi C V^' (k = 1,2). 

Therefore, since L(x) is constant inside the spherical inclusion, we have 

J T(x-y)V l (y)M{ 1} (y)<T(y)dy = v l T l (x-x l )M[ 1) ( T l , (17) 

where x ^ i>i and er^ = (f (y)Vi(y))^ is the random stress field averaged over the volume of the inclusion 
Vi (but not over the ensemble). Hereinafter (y e vj) 

T,(x - xj = (Ui)- 1 J r(x - y)Vi(y)dy, T^x, - x,-) = (T,(y - xj^.. (18) 

The tensor Ty(xj — Xj) has an analytical representation for spherical inclusions of different size in an 
isotropic matrix (see for references Buryachenko, 2007a), the case of ellipsoidal inclusions of different sizes 
and orientations is analyzed by Franciosi and Lebail (2004). 

According to hypothesis HI and to Eshelby's (1961) theorem we get 



where 



<x(x) = BWi(x), vM^aix) = R^(x), x e v h (19) 

B = (I + Q(u i )M^) _1 = const.. R; = ViM^B (20) 
and tensor Q(i>j) is associated with the well-known Eshelby tensor by 

S = I-M (0) Q(w l ). (21) 

Using hypothesis HI and the solution (20) for one inclusion, we get algebraic solution for two 
inclusions V\ and vi subjected to the effective field <Xij(xj): 

2 

Wi(x) = R, 1 ^2 Z..,H,ff,.,;x,i. xetj (i = 1, 2). (22) 
J=l 

Here the matrix Z _1 has the elements (i, j = 1, 2) 

(Z- 1 ),, = My - (1 - -ijR/r^fx, - Xj ). (23) 

For termination of the hierarchy of statistical moment equations in everaged Eq. (9), we will use the 
closing effective field hypothesis: 

H2) Each pair of the inclusions in and Vj is located in an effective field CTjj(x) (22) and 

feW) t = (ff(x)) t ,x6^ (k = i,j). (24) 

After estimating average stresses inside the inclusions (i = 1, 2) 

2 

(<r) l = (c«)- 1 (M«)- 1 ^n«Y ij R j ( CT ), (25) 
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the problem of calculating effective properties of statistically homogeneous and statistically isotropic 
medium becomes trivial and leads (see for details Buryachenko, 2007a) to: 



M* 



where the matrix Y _1 has the following elements (Y -1 )^- = 1,2): 

r 

q=l 

Tij(xi -Xj-)Z|?^(u g ,Xg|;i;i,Xi) -n^T^Xj -x,) <£><:,,■ 



I-R^-^-y T ig (xi -x q )Z^VK,x g |;t; i ,Xi)dx Q 



R; 



(26) 



(27) 



where Z*j = (Z™|) (qr, i,k,l= 1, 2) are the submatrixes (6 x 6) of the matrix Z gi (12 x 12). 

It should be mentioned that the quasi crystalline approximation by Lax (1952) is equivalent to the 
assumption 

Z« = S ZJ I, (28) 
which reduces Eq. (27) (for the aligned homogeneous ellipsoidal particles) to 



(Y-% = %I 



R 



Ty(xj - Xj)yj(u i ,x :7 -|;Ui,x i )-n (: ' ) T i (x i -x,) dx,-, 



(29) 



Due to a widely used second-order probability function Sij(r) (rather than only RDF), we reproduce 
some classic results obtained in the framework of the quasi crystalline approximation by Lax (1952) in 
the form by Willis (1977) (see also Hashin and Shtrikman, 1963; Willis, 1981; Lee, Gillman and Matous, 
2011) 

(<r). = (M^'EY^V), (30) 

3=1 
2 

M* = M (0) + E c^%M[ j) , (31) 



where 



(Y-%=6 ij l- <p®)- l M®Jr(xi - xj)[S y (x< - x,) - cWc^Jdxj, (32) 



To make further progress, the hypothesis of "ellipsoidal symmetry'' 1 for the distribution of inclusions 
attributed to Willis (1977) (see also Khoroshun, 1974, 1978; Ponte Castaneda and Willis, 1995) is widely 
used: 

Hypothesis H3, "ellipsoidal symmetry". The conditional probability density function ip(vj,Xj | 
; Vi, Xj) depends on Xj — Xj only through the combination p = |(ar?.) — (xj — x,)|: 



(p{vj,Xj \;Vi,Xi) = h(p), 



(33) 
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where the matrix (a^-) 1 (which is symmetric in the indexes i and j, a^- = sP^) defines the ellipsoid 
excluded volume = {x : |(a^-) _1 x| 2 < 1}. 

In such a case, the representations of the tensors ~Yij and Yy are simplified 

(Y- 1 ),, = %I-(I-B)^nW, (34) 
(Y _1 )y = y-fe-cWjMfQK ), (35) 

and both pairs of Eqs. (25), (26) and Eqs. (30) (31) are reduced to the same formulae 

(<t>. = B^W+Bc)- 1 ^), (36) 
M* = M (0) + c(cW + Bc^M^B, (37) 

coinciding with the estimations by Mori and Tanaka (1973) method (MTM). 

The destination of the hypothesis H3 is directed towards providing conditions for applying the 
hypothesis HI (see for details Buryachenko, 2010b, 2011). 

The MEFM has a few advantages with respect to the classical one-particle method: 1) Direct binary 
interaction of heterogeneities (see Eq. (22)); 2) Explicit dependence of effective moduli on the partial 
RDF (see Eqs. (25)-(27)); 3) MEFM is not invariant with respect to the size distribution of spherical 
particles (see Eqs. (25)-(27)). In doing so, the classical methods (such as, e.g., Hashin and Shtrikman, 
1963; Hill, 1965; Mori and Tanaka, 1973; and Willis, 1981, see Eqs. (36) and (37)) used the quasi crystalline 
approximation by Lax (1952) are invariant with respect to both the RDF and the size distribution of 
spherical heterogeneities. 



4 Numerical results 



4-1 Unimodal distribution 



The RDF g(r) is well investigated only for identical spherical (3-D and 2D cases, see for references, 
e.g., Buryachenko, 2007a) inclusions with a radius a. Two alternative analytical RDFs of inclusion will 
be examined for 3D case 



g s (r) = H(r-2a) 
g w (r)=H(r-2a) 



1 



2(1 -c 2 



) a 



(38) 
(39) 



where H denotes the Hcaviside step function, r = |x^ — x g | is the distance between the nonintersecting 
fixed inclusions Vi and moving one v q , c is the volume fraction of particles of the radius a. The formula 
(38) describes a well-stirred approximation (differing from the RDF for a Poisson distribution by the 
availability of "excluded volume" with the center where <?(|xi — x 9 |) = 0) while Eq. (39), attributed to 
Willis (1978), takes into account a neighboring order in the distribution of the inclusions (see for details 
I). We will also consider the popular Percus-Yervik (1958) approximation g p ~ Y (r) used for construction 
of Eq. (39): g w (2a) = g p - Y (2a). 

In order to demonstrate the comparison of the available experimental data with the predicting 
capability of the proposed method, we will consider the estimation of the effective elastic moduli L* 
(25). An experimental data set from the work of Smith (1976) used for validation were obtained for 
composite materials consisting of identical glass beads embedded in an epoxy matrix. Young modulus 
and Poisson ratio of the particles and the matrix = 76.0 GPa, i/ 1 ) = 0.23 and E^ = 3.01 GPa, 
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v^°> = 0.394, respectively. As can be seen from Fig. 1 the use of the approach (37) based on the quasi- 
crystalline approximation (27) [also called Mori-Tanaka method (MTM, see Mori and Tanaka, 1973; 
and Benveniste, 1983), curve 4] leads to an underestimate of the normalized effective Young modulus 
E*/E(°> by 1.5 times for c = 0.5 compared with the experimental data. Much better approximations of 
E*/EM are given by the MEFM (23) and (25) with the RDF obtained by the CSM, SCRM, and Willis 
approximation (39) which are distinguished one from another by less than 0.5%. In so doing, the use of the 
well-stirred RDF (38) instead of (39) leads to reduction of predicted values E*/E^ ' on 10%. Therefore, 
in the considered case of the elastic mismatch of constituents, the effective elastic moduli estimated by 
the MEFM (23), (25) are little sensitive with respect the RDF taking into account a neighboring order 
in the particle distribution. In doing so, a reduction of the proposed approach to the classical one (27) 
(which is invariant with respect to any RDFs) leads to significant underestimation of E* /E^ ' . 

Let us now demonstrate an application of the theoretical results by considering an isotropic com- 
posite made of an incompressible isotropic matrix, filled with rigid identical spherical inclusions. This 
example was chosen deliberately because it provides the maximum difference between predictions of rel- 
ative effective shear modulus fi*// 1 i as estimated by the various methods. It should be mentioned that 
in the case being considered, /i*/^ ' is equivalent to relative change of the Newtonian viscosity T]*/r)^ 
of a suspension of identical rigid spheres, see experimental data by Kreger (1972). In Fig. 2 the micromc- 
chanical model MEFM (23) and (25) is analyzed for the effect of choosing different RDFs (sec for details 
I). As can be seen, the effective shear moduli can differ by a factor of two or more depending on the 
chosen RDF. It is interesting that the first peak in Willis (39) approximation is lower than the first peak 
in both the collective rearrangement model (CRM) and shaking collective rearrangement model (SCRM 
simulation, but the RDF (39) decreases slower, and, because of this, the effective shear moduli /i* for the 
RDF (39) larger than fj,* corresponded to the CRM and SCRM. In contrast to the Willis approximation 
(39), PY approximation swiftly decreasing after the first peak leads to the values of //* which are below 
corresponding /i* estimated for both the CRM and SCRM. Just for estimation of the approximation 
quality (171), the curve 2 in Fig. 2 is evaluated for the RDF predicted from a single SCRM's 





0.2 . 0.4 

Concentration of inclusions 



0.0 0.2 . . 0.4 0.6 

Concentration of inclusions 



Fig. 1: E*(c)/EW vs c estimated by the MEFM with 
RDF simulated by CRM (1), SCRM (2), Willis approx- 
imation (3), well-stirred approximation (4), and by the 
MTM (5); o - experimental data by Smith [105] 



Fig. 2. /i*(c)///°) vs c estimated by the MEFM with 
RDF simulated by Willis approximation (1), SCRM's 
approximation (5.5) (2), CRM (3), SCRM (4), PY (5), 
well-stirred approximation (6), and by the MTM (7); o 
- experimental data by Kreger [63] 
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RDF g(r, c) (c = 0.5) according to Eq. (171). As can be seen, the curve 2 provides a reasonably accurate 
approximation of the curve 4 corresponding to the SCRM. The effective moduli fj,* as the functions of c 
(u* ~ c) are very sensitive with respect to the RDF; for example, the use of well-stirred approximation 
(37) leads to reduction of the predicted /x*(c) at c = 0.6 at least by the factor 2 or 3 with respect to the 
RDFs estimated by Eq. (39) and by CRM, and SCRM. Neglecting by the binary interaction of inclusions 
(27) leads to subsequent diminution of /j,*(0.6) in two times. All RDFs taking a neighboring order into 
account lead to infinite values of /x*(c) for large values of c, however the limiting values c ~ 0.61 and 
c ~ 0.62 corresponding the CRM and SCRM, respectively, are better compatible with an experimental 
value c m = 0.637 than the limiting assessments c = 0.58 for the Willis approximation (39). 

4-2 Effective moduli and stress concentrator factors of composites with bimodal packs 

We estimate in Fig. 3 the relative effective shear moduli n* /^°^ of composites with incompressible 
matrix containing rigid spherical inclusions with the RDFs (38) as well as with the RDF simulated by 
the SCRM (I); A 2 i = a^/a^ = 0.313 and &i = c (2) /c (1) = 0.333. 

As can be seen, the values u*/u( ' ~ c estimated by the MEFM for the RDFs corresponding to the 
SCRM's simulation correlates better with experimental data by Chong, Christiansen and Bacr (1971) than 
analogous estimations for the RDF (38). Although the effective behavior of the composite is traditionally 
the main focus of micromechanics, it is also essential to supply insight into the statistical description of 
the local stresses in each phase and at interphase. Estimation of these local fields are extremely useful 
for understanding the evolution of nonlinear phenomena such as plasticity, creep, and damage (see, e.g, 
Ravichandran and Liu, 1995; Tan et al, 2005, 2007). We present Eq. (24) in the form 

<<x) fe = BD fc <<x>, (40) 

where the tensor (A; = 1, 2) (called effective stress concentrator factor) has a simple physical meaning 
of the action of surrounding inclusions on the separate one (<r) fe = Dfc(er). Equation (40) makes it 
possible to estimate the ensemble average of the stress in the vicinity of the inhomogeneities near the 
point x G Vi C v^ k ' with the unit outward normal vector n _L dvi (k = 1, 2) 

(<T-(n)) x = B(n)BD fc ( ( r), (41) 

where the tensor B(n) depends only on the elastic properties of the contacting materials and by the 
direction of the normal n (see for details Buryachenko, 2007a). Thus, the stress concentrators factors 
inside (40) and outside (41) inclusions depend on the linear functions of both the elastic solutions B, 
B(n) for a single inclusion inside infinite matrix and the concentrator factor D& defined intrinsically by 
a micromechanical model. Because of this, we will only analyze a dependence ~ c (k = 1, 2). 

As can be seen in Fig. 4, the components Dunn of the stress concentrator factors (40) for both 
the RDF (38) and SCRM's RDF at the large [k = 1) and small (k = 2) inclusions can be distinguished 
one from another by tens of percent that has a significant practical implementation. Indeed, for periodic 
structures, Zhong and Knaus (2000) numerically demonstrated that for a large difference in particle sizes 
(such as a bimodal distribution) , damage occurs at interfaces between the big particles and matrix, while 
only limited or no damage occurs at interfaces around small ones (that is also experimentally confirmed, 
see Draughn, 1981; Lewandowski, Liu and Hunt, 1989). While these effects of particle interaction and 
size variation are smoothed out in a large ensemble of particles that are exhibited by the effective elastic 
moduli, it is foreseeable that they are an important factor in a failure process such as local debonding 
when the particle size variation leads to a yield-likely response characterized by a nonmonotonic load 
deformation trace. Such a characteristic may, ultimately, be associated with inhomogeneous deformation 
such as occurs in shear bands or other strain concentrations, for example. The effect of increased stress 
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concentrator factors at the big particles for random structure bimodal structures pioneered in the current 
paper is critical for any nonlinear localized phenomena such as onset of yielding, failure initiation, and 
others. Of prime importance is that the mentioned size distribution effect was discovered for linear elastic 
composites with ideal contact between the matrix and particles (analysis of non-ideal contact, see e.g. 
Tan et ai, 2005, 2007, is beyond the scope of the current publication). 




Figure 3: Log[/i*(c)//i' - 1 ] vs c estimated by the 
MEFM with RDFs simulated by the SCRM (curve 1), 
and Eq. (30) (curve 2); o - experimental data by Chong, 
Christiansen and Baer (1971) 



Fig. 4. -Dfc|im vs c estimated for the RDFs simulated 
by the SCRM (curves 1 and 4), and Eq. (26) (curves 2 
and 3) at the big inclusions (1,2) and small ones (3,4) 



5 Conclusion 

5.1 Some critical comments 

It docs not need to be emphasized that popular methods using the second-order probability function 
Sij(r) — 0, 1, . . . ,N) (II) (see for references Buryachenko, 2007a; Lee, Gillman and Matous, 2011) 
are invariant with respect to both £y(r) [at least for statistically isotropic structures, see, e.g., Eqs (36) 
and (37)] and the size distribution of inclusions. Because of this, numerous estimations of Sij(r) (see 
e.g. Kumar, Matous and Geubelle, 2008; Lee et ai, 2009; Lee, Gillman and Matous, 2011; Matous and 
Geubelle, 2006; Liu et ai, 2012) have only a marginal practical use from the perspective of view of 
micromechanics (see Fig. 2 and the comments after Eq. (37)). Moreover, a potential usefulness of Sij(r)'s 
estimations is reduced by a simplified assumption of a two-phase approximation = 0,1) (see. e.g., 
Jackson, Hooks and Buckmaster, 2011) of the 2-point probability density Sij(r) even for multimodal 
distribution (N > 1) of particle sizes. Sij(r) is just one of a few possible statistical descriptors which 
is significantly less sensitive to the structure microtopology than the RDF that provides a smoothness 
of Sij(r) (in opposite to the RDF's estimations which smoothness can be only achieved by additional 
sophisticated efforts, e.g., both the parallelization of simulations and the shaking procedure, see (I)). 
Smoothness Sij(r) and its poor sensitivity in comparison with RDF in analysis of random structures 
provide a possibility for concealment of disadvantages in this simulation. Moreover, coincidence of the 
Sij(r) in both the real sample and the reconstructed structures generated by the protocol dependent 
algorithms (see, e.g., Kumar, Matous and Geubelle, 2008; Lee et ai, 2009) does not assure a coincidence 
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of other more important statistical descriptors (e.g. RDFs). Because of this, comparison of RDFs in both 
the structure generated by protocol independent algorithm (see (I)) and reconstructed structure can be 
considered as a validation of the protocol dependent algorithm used for configuration reconstruction. 
Furthermore, geometrical statistical equivalence of the real and reconstructed packs is necessary but 
not sufficient condition for micromechanical statistical equivalence. In fact, Lee et al. (2009), and Lee, 
Gillman and Matous (2011) implicitly introduced (with an intuitive level of rigor) a new concept of 
Representative Volume Element (RVE) for geometrical statistical equivalence. This concept is not well 
defined and depends on the statistical descriptors used. However, the concept RVE is well known in 
micromechanics for micromechanical statistical equivalence (see for references Buryachcnko, 2007a) which 
is defined not only by Sij (xj — Xj ) (or tp(vj , Xj | ; Vi , Xj ) ) but also by an elastic particle interaction described 
by (p(vj, Xj |; Vi, Xj) rather than by $y(xj — Xj). 

It should be mentioned a growing in importance of advanced experimental techniques (such as X-ray 
tomography and electron microscopy) in analysis of random structures of HSP. In particular, the papers 
referred above (see also for references (I)) demonstrate a large opportunities of statistical descriptions 
of real structures by Sij(xi — Xj). However, a practical using of this statistical descriptors Sij(xi — xj) 
estimated in real HSP is even more important. In this sense it is expected a significant impact on the 
scientific society of the paper Lee, Gillman and Matous (2011) where the authors proposed a direct 
substitution of £y(x, — Xj) estimated in the real packs into the corresponding formulae of analytical 
micromechanics (29) where according to the hypothesis H3 Sij(xi — Xj) should have the ellipsoidal 
symmetry (33). 

A popular explanation of acceptance of this ellipsoidal symmetry hypothesis (33) is that this hy- 
pothesis just simplifies Eq. (32) reducing this equation to Eq. (35) which does not contain the integrals. 
In a similar manner, a destination of the assumption of the ellipsoidal shape of the excluded volume 
in the hypothesis H3 is that this hypothesis just simplifies Eq. (35) by the use of the analytical known 
tensor Q(v®) (expressed through the Eshelby tensor S(vf) (21)) which is exploited instead of a general 
tensor Q(vf)(x) found numerically for an arbitrary shape of v® (see e.g., Subsection 4.7.4 in Buryachenko, 
2007a). However, both mentioned assumptions of the hypothesis H3 have a fundamental conceptual sense 
rather than only an analytical solution of some particular problem. Exploiting the Eshelby tensor con- 
cept in Eq. (4.18) (and in the MEFM) is based on the ellipsoidal shape of the correlation hole w° rather 
than on the inclusion shape Wj. An abandonment of cither the assumption of the v^s ellipsoidal shape or 
ellipsoidal symmetry hypothesis (32) necessarily leads to the inhomogencity of the effective field <Xi(y) 
(y £ v i) acting on the inclusion Vi that is prohibited for the classic background of micromechanics based 
on Eq. (9) (see for details Subsection 5.3 and Buryachenko, 2010a, 2010b, 2010c, 2012). 

Thus, correct performing of micromechanical analysis in Section 3 uniquely determines the form of 
statistical descriptors either (p(vj,Xj\;Vi,Xi) or Sij(xi — Xj) estimated from the real packs. At first, it is 
necessary to check that the real pack has a statistically anisotropic structure rather than a functional 
graded configuration (as in Lee, Gillman and Matous, 2011) requiring a different level of micromechanical 
modeling (see for details Buryachenko, 2007a). For statistically anisotropic packs, the mentioned descrip- 
tors must be approximated by the appropriate functions with ellipsoidal symmetry (33), otherwise the 
micromechanical models corresponding to classic background of micromechanics (9) can not be used (al- 
though the new background of micromechanics considered in Subsection 5.3 makes it possible to exploit 
any descriptors). In any way, ip(vj,Xj\;vi,Xi) is preferred over S^x, — Xj) due to better sensitivity to a 
concrete microstructure (either statistically isotropic, statistically anisotropic, or functionally graded) as 
well as to exploiting precisely <p(vj, Xj\; Vi, Xj) (rather than Sij(xi — Xj)) in the advanced micromechanical 
methods (see Fig. 2 and Buryachenko, 2007a, 2011). 

5.2 Some prospective problems 



V.A. Buryachenko 



13 



The slightly modified version of the MEFM proposed provides the calculation with reasonable accu- 
racy of the effective elastic moduli and statistical average of stresses in the components for the composites 
with bimodal distribution of spherical particles. A straightforward generalization of this model is possible 
to the multimodal particle distribution. It also allows the model to simultaneously account for a plethora 
of constituents in the system - for example, metal powder, short fibers, and voids can be included. 

However, the most promising result of the approach proposed is sensitivity of estimations of both 
the effective elastic moduli and stress concentrator factors on the size ratios of particles that is espe- 
cially critical for the HSP with extremely high volume fractions (> 0.9). A similar effect is expected to 
be detected for other linear problems such as e.g. thermoconductivity, thermoelasticity, viscoelasticity 
(Banerjee and Adams, 2004) which are concerned with a number of propellant-related issues, including 
the burning rate and the thermomechanical loading of the HSPs. Intensively investigated micromechani- 
cal problem with imperfect interface between the matrix and particles are solved by substitution into the 
known cohesive zone model the statistical averages of the local stresses at the oxidizer particles estimated 
by the Mori-Tanaka method (see for references Tan et al., 2005, 2007). These estimations (depending 
only on the volume fraction c rather than on the RDFs gij and the size ratios A,j of particles) can be 
significantly refined by the proposed version of the MEFM. 

Evolution of nonlinear phenomena such as viscoplasticity, creep, and damage are a new challenge 
for designers of rockets. When one tries to estimate the equivalent stress in the strength theories as well 
as in nonlinear creep theory, or when the yield function in plasticity theory is considered, squares of the 
first invariant or the second invariant of the deviator of local stresses are frequently used. For unimodal 
distribution of particles, estimations of the required second statistical moments of local stresses can be 
performed by both the perturbation method and the method of integral equations (see for details Chapter 
13 in Buryachenko, 2007a). Generalization of the mentioned methods to the case of the multimodal 
distribution of particles (which is inherent in the HSP) is straightforward. It opens a new avenue of 
questions on some specific nonlinear problems such as ignition, hot spot analysis, and detonation. So 
Field et al. (1992), Massoni et al., 1999; Tan et al. (2005), Dienes, Zuo, and Kershner (2006), Grinfeld 
(2009) indicated the following mechanisms of energy concentrations in locations called "hot-spot" causing 
localized melting, ignition, and fast burn of the reactive material: 1) Adiabatic compression of cavity 
gases; 2) Heating of solid adjacent to collapsing cavity; 3) Viscous heating of binder between grains; 
4) Friction between impacting surfaces; 5) Localized adiabatic shear; 6) Heating by frictional sliding 
between the faces of a closed cracks and debondings; 7) Debonding of particles and matrix. All of these 
problems were partially solved for HSPs and explosives by simplified micromechanical models such as, 
e.g., the dilute approximation ((<?},. = (<x)), the effective medium approximation considering each defect 
inside the effective medium (i/ ' = L*), and the mixture role of damage mechanics. It can also be 
mentioned that substitution of estimations of both the first and second statistical moments of local stresses 
performed by the proposed approach (see Sections 3 and 4) into the mentioned nonlinear problems is also 
straightforward. 

However, more detailed consideration of particular problems mentioned is beyond the scope of the 
current paper and will be considered in subsequent publications. More rich in content is a discussion of the 
main hypotheses as well as the limitations of the proposed estimations and their possible generalizations 
for solutions of some prospective problems. 

5.3 Generalization of some hypotheses 

The assumption HI of homogeneity of <r(x), (x G is a classical hypothesis of micromcchanics (see 
for the earliest references Lax, 1952) and was accepted in order to make it easier to solve the problems 
for both one (19) and two (22) particles. Buryachenko (2007a) demonstrates that the MEFM includes 
plurality of popular methods based on the hypothesis HI. The accuracy of this hypothesis was estimated 
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for two spherical inclusions in the infinite matrix as well as for the periodic structure composites (see for 
references Buryachenko, 2007a). 

However, the hypothesis HI) is merely a zero-order approximation of binary interacting inclusions 
that results in a significant shortcoming of the MEFM. This substantial obstacle can be overcame in the 
upgraded version of the MEFM proposed by Buryachenko (2007b) in light of the generalized schemes 
based on the numerical solution (some sort of the building blocks) of the problem for both one and two 
inclusions in the infinite media, subjected to the effective field evaluated from forthcoming self-consistent 
estimations. It is possible to get a numerical solution in the form of a table in which each line among 
a few hundred lines corresponds to an accurate solution for two inclusions with concrete coordinates 
in an infinite medium subjected to the unit loading. No restrictions are imposed on the microtopology 
of the microstructure and the shape of inclusions as well as on the inhomogeneity of stress field inside 
the inclusions. However, the main computational advantage of the generalized MEFM lies in the fact 
that such fundamental notions of micromechanics as a Green function and Eshelby tensor are not used, 
and we can analyze any anisotropy of constituents (including the matrix) as well as any shape and any 
composite structure of inclusions. The known numerical methods such as FEA, VIE, BIE, and the multi, 
multipole expansion method (see e.g., Hatch and Davis, 2006; Buryachenko, Kushch, and Roy, 2007), and 
complex potential method (Buryachenko and Kushch, 2006) which can be used for construction of the 
building blocks mentioned, have a series of advantages and disadvantages. It is crucial for the analyst to 
be aware of their range of applications. In particular, for the high volume fraction of particles c, phase 
arrangements will involve frequent direct contacts between neighboring particles and may even lead to 
a local stress singularities between the contacting rigid particles. However, these strong inhomogeneity 
of stress distributions can only be detected and analyzed by advanced accurate numerical methods (see, 
e.g., Hatch and Davis, 2006; Buryachenko and Kushch, 2006) while the analytical method (22) used in 
the current paper takes into account only the stresses averaged over the volumes of neighboring particles. 
The accuracy of the method (22) estimated in Fig. 4.4 in Buryachenko (2007a) eliminates the possible 
difficulties which may arise when direct contact between neighboring particles occurs frequently. 

Fundamentally more general approach is based on a new general integral equation of micromechanics 
proposed by Buryachenko (2010a, 2010b) 



Equation (34) was obtained without any auxiliary assumptions such as, e.g., the version of the HI: 
(r(x— y)Tj(y))(y) = r(x — y)(r;)(y) (hypothesis Hlb, p. 253 in Buryachenko, 2007a) implicitly exploited 
in the known centering methods and reducing Eq. (42) for statistically homogeneous media subjected to 
the homogeneous boundary conditions to the known Eq. (9) which goes back to Lord Rayleigh (1892). 
Buryachenko (2010a, 2010b) demonstrated that Eq. (9), erroneously recognized as an exact one after the 
proofs by Shermergor (1977), is correct only after the additional asymptotic assumption Hlb. What seems 
to be only a formal trick is in reality a new background of micromechanics allowing us to abandon the 
central concept of classical micromechanics, such as effective field hypothesis HI. Then the hypothesis 
H2 taking into account the binary interaction of inclusions leads to the inhomogeneity of both the 
statistical average stresses inside the inclusions and the effective fields (violation of the effective field 
hypothesis HI) even for statistically homogeneous composites subjected to the homogeneous remote 
loading and containing homogeneous ellipsoidal inclusions. It is expected significant improvement of 
accuracy of statistical averages of local stress concentrator factors (with a possible change of the sign of 
predicted local stresses) in the approach by Buryachenko (2011) (only applied for 2D case) with respect 
to a classical one presented in this paper while their averaged values (cr) i (and, therefore, the effective 
moduli L*) arc only slightly sensitive to the random fiber arrangement described by the RDF. 




(42) 
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Here a question of verification of the approach proposed appears. Unfortunately, validating the nu- 
merical results expected is currently beyond the capabilities of standard experimental techniques based 
on the measurement of effective elastic moduli L* because the difference of L* estimated by the new 
Buryachenko (2011) and Buryachenko (2007a) approaches is expected negligible (at least for the moder- 
ate volume fraction c). However, the advantage of the method proposed can be seen in the fact that the 
emphasis in materials modeling shifted from explaining experimentally observed properties of materials 
(such as L*) to predicting not-yet-measured properties (such as the local stresses). Therefore, there is 
an increasing need to provide very detailed experimental information about the local stress concentrator 
factors of numerically investigated materials at the microscale. New volumetric digital image correlation 
(VDIC) relying on the X-ray tomographic imaging of naturally occurring material texture within samples 
requires the ability to reproduce the position of image points with "high precision" because the determi- 
nation of relative deformations requires two images for comparison purposes to extract strain measures. 
Full-field measurements of strains by VDIC (and by other imaging tools) are potentially well suited to 
analyze the specific mechanical properties of composite materials (see, e.g., Grcdiac, 2004; Sutton, Ortcu 
and Schreier, 2009). However, the mentioned opportunities of experimental methods are not yet realized 
at the level required for verification of the new approach mentioned above. 

Of course, all particular problems mentioned in Subsections 5.1 and 5.2 can be solved in the frame- 
work of the new background of micromechanics indicated in the current subsection 5.3 (see also Bury- 
achenko, 2011, 2012). 
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